Carga de librerías

if (!(require(DT)))
  install.packages("DT")
library(DT)

if (!require(readxl))
  install.packages("readxl")
library(readxl)

# Rfast para matriz de varianza combinada
if (!requireNamespace("Rfast"))
  install.packages("Rfast")
library(Rfast)

if (!requireNamespace("ggplot2"))
  install.packages("ggplot2")
library(ggplot2)

if (!requireNamespace("plotly"))
  install.packages("plotly")
library(plotly)

if (!requireNamespace("ggrepel"))
  install.packages("ggrepel")
library(ggrepel)

# GridFCM
library(devtools)
if (!requireNamespace("GridFCM.practicum"))
  install_github("asanfe/GridFCM.practicum", quietly = TRUE)
library(GridFCM.practicum)

# Viridislilte
if (!requireNamespace("viridisLite"))
  install.packages("viridisLite")
library(viridisLite)

# Test para normalidad multivariante
if (!requireNamespace("MVN"))
  install.packages("MVN")
library(MVN)

if (!requireNamespace("ggpattern", quietly = TRUE))
  install.packages("ggpattern")
library(ggpattern)

if (!requireNamespace("factoextra", quietly = TRUE))
  install.packages("factoextra")
library(factoextra)

if (!requireNamespace("cluster", quietly = TRUE))
  install.packages("cluster")
library(cluster)

if (!requireNamespace("RColorBrewer", quietly = TRUE))
  install.packages("RColorBrewer")
library(RColorBrewer)

Importación de WIMP

path <- 'Wimp_Ejemplo.xlsx'
opr <- TRUE

wimp <- GridFCM.practicum::importwimp(path = path, opr = opr, sheet = 1)
bertin(wimp$openrepgrid, colors = c("palegreen", "darkgreen"))

Digrafo con GridFCM.practicum

# Digraph
GridFCM.practicum::digraph(wimp, layout = "rtcircle")

E/S de los constructos. Método Simple

c.io.test <- GridFCM.practicum::degree_index(wimp, method = "simple")
c.io.test <- c.io.test[, c(2,1,3)]
c.io.test <- cbind(c.io.test, Diff = (c.io.test[, 2] - c.io.test[, 1]))
c.io.test.r <- round(c.io.test, 3)
DT::datatable(c.io.test.r)

E/S de los constructos. Método wnorm

c.io.test <- GridFCM.practicum::degree_index(wimp, method = "wnorm")
c.io.test <- c.io.test[, c(2,1,3)]
c.io.test <- cbind(c.io.test, Diff = (c.io.test[, 2] - c.io.test[, 1]))
c.io.test.r <- round(c.io.test, 3)
DT::datatable(c.io.test.r)

Ejemplo de salida de P-H index

test.wphm <- GridFCM.practicum::ph_index(wimp = wimp, method = "weight", std = FALSE)
DT::datatable(test.wphm)

Distancia de Mahalanobis y distribución de datos

Test de Mardia para análisis multivariante

Llevamos a cabo previamente un test de Mardia para constrastar la normalidad multivariante de los datos, a fin de determinar la pertinencia del punto de corte basado en adecuación a distribución Chi-cuadrado de distancia de Mahalanobis

# Test de Mardia

test.result <- mvn(data = test.wphm, mvnTest = "mardia")

print(test.result)
## $multivariateNormality
##              Test          Statistic           p value Result
## 1 Mardia Skewness   5.98339909149846 0.200391474273193    YES
## 2 Mardia Kurtosis -0.360511421430716 0.718464716991934    YES
## 3             MVN               <NA>              <NA>    YES
## 
## $univariateNormality
##               Test  Variable Statistic   p value Normality
## 1 Anderson-Darling     p        0.4492    0.2505    YES   
## 2 Anderson-Darling     h        0.8801    0.0198    NO    
## 
## $Descriptives
##    n         Mean  Std.Dev     Median       Min      Max       25th      75th
## p 21 4.806081e+00 2.341929  4.2544258  1.308148 9.981991  3.4530381 5.9279118
## h 21 2.641323e-17 1.694968 -0.2946278 -3.335187 3.346972 -0.7896026 0.4596194
##         Skew   Kurtosis
## p 0.70632887 -0.2634039
## h 0.08923869 -0.2668683

Test de resultado de la función

test.wmahalanobis <- GridFCM.practicum::mahalanobis_index(wimp = wimp, method = "weight", std = FALSE, sign.level = 0.2)
DT::datatable(test.wmahalanobis)

Gráfica de barras de distancias de Mahalanobis y punto de corte

# Distancia de Mahalanobis
test.bp.wmahalanobis <- GridFCM.practicum::mahalanobis_index(wimp = wimp, method = "weight", std = FALSE)
test.wmahalanobis.df <- as.data.frame(test.bp.wmahalanobis)


#test.wmahalanobis.df$constructo <- rownames(test.wmahalanobis)
test.wmahalanobis.df$constructo <- wimp$constructs$right.poles

# Valoración del ideal
test.wmahalanobis.df$idealdirect <- wimp$ideal$direct

# Columna para identificar constructos dilemáticos
test.wmahalanobis.df$fill.color <- ifelse(test.wmahalanobis.df$idealdirect == 4, "yellow2", "honeydew")

# Ordenamos las barras en orden decreciente
test.wmahalanobis.df <- test.wmahalanobis.df %>%
  arrange(desc(m.dist))

# Convertimos 'constructo' en un factor con los niveles en el orden deseado
test.wmahalanobis.df$constructo <- factor(test.wmahalanobis.df$constructo, levels = test.wmahalanobis.df$constructo)

# Punto de corte distribución Chi-Cuadrado
sign.level <- 0.2
df <- ncol(test.wphm)
chi.square.cutoff <- qchisq(1 - sign.level, df)
#media_m_dist <- mean(test.wmahalanobis.df$m.dist)

Constructos supraordenados, resaltando dilemáticos

# Crear el histograma de constructos supraordenados

# Filtramos por los constructos donde el valor de 'h' es mayor que cero
test.wmahalanobis.df.sup <- test.wmahalanobis.df %>%
  filter(h > 0)

bar_plot <- ggplot(test.wmahalanobis.df.sup, aes(x = constructo, y = m.dist, fill = fill.color)) +
  geom_bar(stat = "identity", color = "black", linewidth = 0.25) +
  geom_hline(yintercept = chi.square.cutoff, linetype = "dashed", color = "darkgreen", linewidth = 1) +
  scale_fill_identity() + # Usa los colores asignados directamente
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "none" # Ocultar leyenda para fill
  ) +
  labs(x = "Constructos con h > 0", y = "Distancia de Mahalanobis", title = "Constructos con h>0 por distancia de Mahalanobis")

# Mostramos el gráfico
print(bar_plot)

Constructos subordinados, resaltando dilemáticos

# Crear el histograma de constructos subordinados

# Filtramos por los constructos donde el valor de 'h' es menor que cero
test.wmahalanobis.df.sub <- test.wmahalanobis.df %>%
  filter(h < 0)

bar_plot <- ggplot(test.wmahalanobis.df.sub, aes(x = constructo, y = m.dist, fill = fill.color)) +
  geom_bar(stat = "identity", color = "black", linewidth = 0.25) +
  geom_hline(yintercept = chi.square.cutoff, linetype = "dashed", color = "darkgreen", linewidth = 1) +
  scale_fill_identity() + # Usa los colores asignados directamente
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "none" # Ocultar leyenda para fill
  ) +
  labs(x = "Constructos con h < 0", y = "Distancia de Mahalanobis", title = "Constructos con h<0 por distancia de Mahalanobis")

# Mostramos el gráfico
print(bar_plot)

Ejemplo de representación en espacio P-H

Función de representación ggplot

graph_ph_gg <- function(ph.mat){
  
  # Convertimos matriz en dataframe
  ph.mat.df <- as.data.frame(ph.mat)
  
  # Añadimos los nombres de los constructos como una nueva columna en el dataframe
  ph.mat.df$constructo <- rownames(ph.mat)
  
  # Paleta de colores en escala de verdes
  paleta.verdes <- viridis(n = length(ph.mat.df$constructo), option = "viridis")
  
  # Gráfica de dispersión con tema de fondo blanco
  wimp.plot <- ggplot(ph.mat.df, aes(x = p, y = h)) +
    geom_point(aes(fill = constructo), shape = 21, size = 3, color = "black") +
    geom_text_repel(aes(label = constructo), size = 3.5, fontface = "plain", box.padding = 0.5) +
    scale_fill_manual(values = paleta.verdes) +
    geom_abline(slope = -1, intercept = 0, linetype = "dashed", color = "darkgreen") +
    geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "darkgreen") +
    labs(
      x = "P - Presencialidad (frecuencia del constructo)",
      y = "H - Jerarquía (influencia del constructo)",
      title = "Gráfica de Dispersión de Constructos en el Espacio P - H"
    ) +
    theme_bw() + 
    theme(
      panel.grid.major = element_line(linewidth = 0.5, linetype = 'solid', colour = "lightgrey"),
      panel.grid.minor = element_blank(), 
      legend.position = "none"  
    ) +
    xlim(c(0, max(ph.mat.df$p))) 
    #ylim(c(-0.5, 0.5))

  return(wimp.plot)  
}

Representación de espacio PH

Sin estandarización - ggplot

test.wphm.se <- GridFCM.practicum::ph_index(wimp = wimp, method = "weight", std = FALSE)

wp1 <- graph_ph_gg(test.wphm.se)
print(wp1)

Sin estandarización - plotly sin marcar área no viable ni constructos centrales

wp1p <- graph_ph(test.wphm.se, mark.nva = FALSE, mark.cnt = FALSE)

wp1p

Espacio PH con coloreado de área no útil y marcado de outliers. Función de representación

Sin estandarización

# Matriz PH con distancia de Mahalanobis
test.bp.wmahalanobis.se <- GridFCM.practicum::mahalanobis_index(wimp = wimp, method = "weight", std = FALSE, sign.level = 0.2)

wp3 <- graph_ph(test.bp.wmahalanobis.se, mark.nva = TRUE, mark.cnt = TRUE)

wp3

Sin estandarización - plotly sin marcar area de coordenadas no viables. Sin puntos

wp3.sp <- graph_ph(test.bp.wmahalanobis.se, mark.nva = TRUE, mark.cnt = TRUE, show.points = FALSE)

wp3.sp
## Warning: No trace type specified and no positional attributes specified
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Sin estandarización - plotly sin marcar area de coordenadas no viables

wp1mhp <- GridFCM.practicum::graph_ph(test.bp.wmahalanobis.se, mark.nva = FALSE, mark.cnt = TRUE)

wp1mhp

Primer eigenvectorsobre matriz de implicaciones

eigen_index <- function(wimp){

  # Matriz de adyacencia
  adj.matrix <- wimp$scores$implications

  # Vectores y valores propios
  results <- eigen(adj.matrix)

  # Extrae el primer vector propio, asociado a la varianza máxima de los datos
  eigenvector.principal <- results$vectors[,1]

  # Creamos un marco de datos con los nombres de los constructos y las componentes del primer eigenvector
  df.centrality <- data.frame(
    constructs = wimp$constructs$constructs,
    leftpoles = wimp$constructs$left.poles,
    rightpoles = wimp$constructs$right.poles,
    firsteigenvector = eigenvector.principal
  )

  return(df.centrality)
}
cent.evalues.df <- eigen_index(wimp)
DT::datatable(cent.evalues.df)

Análisis de componentes principales

# Foco del PCA
adj.matrix <- wimp$scores$implications

# Análisis de componentes principales
pca.result <- prcomp(adj.matrix, center = TRUE, scale = TRUE)

# Extraer los dos primeros componentes principales
pca.comp <- as.data.frame(pca.result$x[, 1:2])
pca.comp$constructs <- wimp$constructs$constructs

# Crea la gráfica de dispersión
pca.plot <- plot_ly(data = pca.comp, x = ~PC1, y = ~PC2, type = 'scatter', mode = 'markers',
                    hoverinfo = 'text+x+y',
                    marker = list(size = 10, opacity = 0.8)) %>%
            layout(title = 'PCA de matriz de adyacencia',
                   xaxis = list(title = 'PCA 1'),
                   yaxis = list(title = 'PCA 2'),
                   hovermode = 'closest',
                   plot_bgcolor = "white",
                   font = list(family = "Arial"),
                   showlegend = FALSE) %>%
            # Add annotations for each point
            add_annotations(data = pca.comp, x = ~PC1, y = ~PC2, text = ~constructs,
                            showarrow = FALSE, xanchor = 'center', yanchor = 'bottom', font = list(size = 12))

# Muestra la gráfica
pca.plot

Análisis por conglomerados

Determinación de número óptimo de conglomerados

# Coordenadas de los constructos en espacio P-H
test.dist <- GridFCM.practicum::ph_index(wimp = wimp, method = "weight", std = FALSE)
rownames(test.dist) <- as.character(wimp$constructs$right.poles)

# Matriz de distancias entre los constructos
distancias <- dist(test.dist, method = "euclidean")

# Cáclulo de número de conglomerados por coeficiente de la silueta-----
# Vector para guardar los promerios de los coeficientes de silueta
sil.widths <- numeric()

# Número máximo de clusters a evaluar. Acotamos, como máximo, al número de constructos disponible menos 1 
max.clusters <- length(wimp$constructs$constructs) - 1

# Calculamos PAM y el coeficiente de silueta para diferentes números de clusters
for(j in 2:max.clusters) {
  pam.stp <- pam(distancias, j)
  sil.stp <- silhouette(pam.stp)
  sil.widths[j] <- mean(sil.stp[, "sil_width"])
}

# Identificar el número de clusters que maximiza el coeficiente de silueta
k <- which.max(sil.widths)
k
## [1] 2

Tenemos un número máximo de 2 conglomerados en nuestros datos, obteniendo en ese caso un valor medio de anchura de silueta de 0.5660763

Representación de números de conglomerados óptimo

Adecuación de cohesión y separación de cada punto según pertenezca a distintos conglomerados:

# Lista que albergará las distintas gráficas de silueta
lista.graf.sil <- list()

# Preparamos una lista con diversas representaciones gráficas de siluetas (de 2 a 10 clústeres)
for(j in 2:13){
  it.pam <- cluster::pam(distancias, j, diss = TRUE)
  p <- factoextra::fviz_silhouette(it.pam, label = FALSE, print.summary = FALSE)
  lista.graf.sil[[j-1]] <- p
}

# Organizar los gráficos en una matriz de 4x3, y los presentamos
gridExtra::grid.arrange(grobs = lista.graf.sil, ncol = 3, nrow = 4)

Dendrograma

# Paleta de colores en escala de verdes
#paleta.verdes <- viridis(n = k, option = "viridis")
#paleta.verdes <- rev(brewer.pal(n = k, name = "Greens"))
#paleta.verdes <- sapply(paleta.verdes, function(color) RColorBrewer::adjustcolor(color, brightness = 0.5))
paleta.verdes <- c("#003300","#008000", "#3CB371")

#Dendrograma
modelo.hclust <- hcut(distancias, k = k, stand = TRUE)
fviz_dend(modelo.hclust, rect = TRUE, cex = 0.7,
          k_colors = paleta.verdes,horiz = TRUE,
          main = 'Dendrograma de constructos por distancia en P-H')
## Registered S3 method overwritten by 'dendextend':
##   method       from   
##   text.pvclust pvclust
## Warning in get_col(col, k): Length of color vector was longer than the number
## of clusters - first k elements are used
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.